library(dplyr)
library(stringr)
library(ggplot2)
library(psych)
library(ggthemes)
library(estimatr)

jordan = read.csv("jordan_svy.csv", stringsAsFactors = F)

#load pca function and coefplot function
source("pca_function.R")
source("coefplot_function.R")


# Data preparation --------------------------------------------------------
jordan$treatment = factor(jordan$treatment, levels = c("0", "1", "2", "3"), labels = c("Control", "Values", "Need", "Joy"))

jordan = jordan %>% 
mutate(out_border = 8 - out_border, #For these outcomes, I'm changing it so that scoring higher would mean having better attitudes toward refugees
       out_quarantine = 8 - out_quarantine,
       out_sendback = 8 - out_sendback,
       out_volunteer1 = 2 - out_volunteer1, #Changing it so that 1 means agreed to volunteer and 0 did not agree
       out_volunteer2 = 2 - out_volunteer2) %>% 
  mutate(contact = ifelse(interact == 2, 0, #frequency of contact with Syrians
                          ifelse(interact_freq %in% c(6, 7), NA, 6 - interact_freq)), #Reverses coding of interaction frequency.
         income_low = ifelse(income < 4, 1, 0),
         unemployed = ifelse(econ_sit == 3, 1, 0),
         religious = ifelse(gender == 1 & mosque == 1, 1, #Measuring it just for males since females are very unlikely to attend mosque 
                            ifelse(gender == 1 & mosque > 1, 0, NA)), 
         westbank = ifelse(origin2=="west", 1,
                           ifelse(origin2=="east",0,NA)),
         low_educ = ifelse(educ == 1, 1,
                           ifelse(educ > 1, 0, NA)),
         low_knowledge = ifelse(knowledge < 3, 1,
                                ifelse(knowledge >= 3, 0, NA)),
         high_contact = ifelse(contact > 2, 1, ifelse(contact <= 2, 0, NA)),
         high_density = ifelse(refugees_num > quantile(refugees_num, probs = .8), 1,
                               ifelse(refugees_num < quantile(refugees_num, probs = .8), 0, NA))
  )

# pca ---------------------------------------------------------------------
all = c("out_host", "out_thermo", "perceive_rlns", "out_housing", "out_econ", "out_crime", "out_impact", "out_services", "out_terror", "out_educ", "out_culture", "out_agri", "out_image", "out_quarantine", "out_workpermits", "out_sendback", "out_border", "out_volunteer1", "out_volunteer2")
attitudes = c("out_host", "out_thermo", "perceive_rlns")
impact = c("out_housing", "out_econ", "out_crime", "out_impact", "out_services", "out_terror", "out_educ", "out_culture", "out_agri", "out_image")
policy = c("out_quarantine", "out_workpermits", "out_sendback", "out_border")
volunteer = c("out_volunteer1", "out_volunteer2")

jordan2 = jordan

#making sure behavior has no NAs

jordan_all = jordan2[complete.cases(jordan2[, all]),]
jordan_attitudes = jordan2[complete.cases(jordan2[, attitudes]),]
jordan_impact = jordan2[complete.cases(jordan2[, impact]),]
jordan_policy = jordan2[complete.cases(jordan2[, policy]),]
jordan_volunteer = jordan2[complete.cases(jordan2[, volunteer]),]

jordan_all = jordan_all %>% as.data.frame()
jordan_attitudes = jordan_attitudes %>% as.data.frame()
jordan_impact = jordan_impact %>% as.data.frame()
jordan_policy = jordan_policy %>% as.data.frame()
jordan_volunteer = jordan_volunteer %>% as.data.frame()



#get PCs
jordan_all$pca_all = PCA_extract(jordan_all[, all], cor_type = "pear")[[2]] %>% as.vector()

jordan_attitudes$pca_attitudes = PCA_extract(jordan_attitudes[, attitudes], cor_type = "mix")[[2]] %>% as.vector()

jordan_impact$pca_impact = PCA_extract(jordan_impact[, impact], cor_type = "poly")[[2]] %>% as.vector()

jordan_policy$pca_policy = PCA_extract(jordan_policy[, policy], cor_type = "poly")[[2]] %>% as.vector()

jordan_volunteer$pca_volunteer = PCA_extract(jordan_volunteer[, volunteer], cor_type = "poly")[[2]] %>% as.vector()

# Coef plots on PCAs ------------------------------------------------------
names = c("Age (years)", "Education (1-6)", "Income (1-12)", "West Bank (binary)", "Self-Reported Contact (0-5)", "No. of Refugees in District (log)", "Location: north (binary)", "Location: south (binary)", 
          "Religiosity (0-7)", "Cultural Threat (PC)", "Knowledge (1-5)", "Female", "Household (1-12)")

model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all)
model2 = lm_robust(pca_attitudes ~ treatment, weights = weights, data = jordan_attitudes)
model3 = lm_robust(pca_policy ~ treatment, weights = weights, data = jordan_policy)
model4 = lm_robust(pca_impact ~ treatment, weights = weights, data = jordan_impact)
model5 = lm_robust(pca_volunteer ~ treatment, weights = weights, data = jordan_volunteer)

model1 = tidy(model1)
model2 = tidy(model2)
model3 = tidy(model3)
model4 = tidy(model4)
model5 = tidy(model5)

model1 %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ outcome)
ggsave("figures/figureC1.pdf", width = 4, height = 4)

model_final = bind_rows(model2, model3, model4, model5)

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(outcome = recode(outcome, "pca_attitudes" = "Attitudinal Outcomes",
                          "pca_policy" = "Policy Outcomes", "pca_impact" = "Impact on Jordan",
                          "pca_volunteer" = "Volunteer Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ outcome)
ggsave("figures/figureC2.pdf", width = 6, height = 6)


# Competition ---------------------------------------------------
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$income_low==1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$income_low==0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Below Median Income", "2" = "Median Income or Above")) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC3.pdf", width = 5, height = 3)


model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$unemployed==1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$unemployed==0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Unemployed", "2" = "Other Job Status"),
         model = factor(model, levels = c("Unemployed", "Other Job Status"))) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model, scales = "free")
ggsave("figures/figureC4.pdf", width = 5, height = 3)

# Religiosity -------------------------------------------------------------
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$religious == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$religious == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% filter(term == "treatmentValues") %>%  mutate(term = str_remove_all(term, "treatment")) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Attend Mosque Daily", "2" = "Attend Mosque \nLess Often"),
         model = factor(model, levels = c("Attend Mosque Daily", "Attend Mosque \nLess Often"))) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = model, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Religiosity", subtitle = "", y = "Estimate")
ggsave("figures/figureC5.pdf", width = 5, height = 3)


# shared experience -------------------------------------------------------
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$westbank == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$westbank == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% filter(term == "treatmentNeed") %>%  mutate(term = str_remove_all(term, "treatment")) %>% 
  filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "West Bank Origin", "2" = "East Bank Origin"),
         model = factor(model, levels = c("West Bank Origin", "East Bank Origin"))) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = model, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Origin", subtitle = "", y = "Effect of Need Treatment")
ggsave("figures/figureC6.pdf", width = 5, height = 3)



# Education/Sophistication ---------------------------------------------------------------
#Education
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$low_educ == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$low_educ == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Below Median Education", "2" = "Median Education or Above")) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC7.pdf", width = 5, height = 3)


#Political Knowledge
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$low_knowledge == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$low_knowledge == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Below Median \nPolitical Knowledge", "2" = "Median Political \nKnowledge or Above")) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC8.pdf", width = 5, height = 3)


# Contact -----------------------------------------------------------------
#Self reported
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$high_contact == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$high_contact == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Above Median Contact", "2" = "Below Median Contact")) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC9.pdf", width = 5, height = 3)


#High density
model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$high_contact == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$high_density == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "High Refugee Density", "2" = "Low Refugee Density")) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC10.pdf", width = 5, height = 3)



# Young, urban, educated --------------------------------------------------
jordan_all = jordan_all %>% mutate(online_sample = ifelse(age < 40 & location == "center" & low_educ == 0, 1, 0))

model1 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$online_sample == 1,]) %>% tidy
model2 = lm_robust(pca_all ~ treatment, weights = weights, data = jordan_all[jordan_all$online_sample == 0,]) %>% tidy

model_final = bind_rows(model1, model2, .id = "model")

model_final %>% mutate(term = str_remove_all(term, "treatment")) %>% filter(term!="(Intercept)") %>% 
  mutate(model = recode(model, "1" = "Young, Urban, and \nEducated Respondents", "2" = "Other Respondents"),
         model = factor(model, levels = c("Young, Urban, and \nEducated Respondents", "Other Respondents"))) %>% 
  mutate(outcome = recode(outcome, "pca_all" = "All Outcomes")) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange()  + 
  theme_few() + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  labs(x = "Treatment", subtitle = "", y = "Estimate") + 
  facet_wrap(~ model)
ggsave("figures/figureC11.pdf", width = 5, height = 3)
